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Abstract 

The position-dependent exact-exchange energy per particle e x (z) (defined as the interaction 
between a given electron at z and its exact-exchange hole) at metal surfaces is investigated, by using 
either jellium slabs or the semi-infinite (SI) jellium model. For jellium slabs, we prove analytically 
and numerically that in the vacuum region far away from the surface e^ lah (z — > oo) — > — e 2 /2z, 
independent of the bulk electron density, which is exactly half the corresponding exact-exchange 
potential V x (z -> oo) -> - e 2 /z [Phys. Rev. Lett. 97, 026802 (2006)] of density-functional theory, 
as occurs in the case of finite systems. The fitting of e^ lab (z) to a physically motivated image- 
like expression is feasible, but the resulting location of the image plane shows strong finite-size 
oscillations every time a slab discrete energy level becomes occupied. For a semi-infinite jellium, 
the asymptotic behavior of e^ l (z) is somehow different. As in the case of jellium slabs £ x l (z — > oo) 
has an image- like behavior of the form oc — e 2 /z, but now with a density-dependent coefficient that 
in general differs from the slab universal coefficient 1/2. Our numerical estimates for this coefficient 
agree with two previous analytical estimates for the same. For an arbitrary finite thickness of a 
jellium slab, we find that the asymptotic limits of e^} ah (z) and £ x l {z) only coincide in the low-density 
limit (r s — > oo), where the density-dependent coefficient of the semi-infinite jellium approaches the 
slab universal coefficient 1/2. 
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I. INTRODUCTION 



The jellium model of a metal surface, introduced by Bardeen in 1936, 1 is the simplest 
model which reproduces qualitatively, and sometimes quantitatively, the physical properties 
of real metal surfaces.- While in his work Bardeen applied an approximated Hartree-Fock 
(HF) theory for the study of the electronic structure, since the seminal work of Lang and 
Kohn (LK)2 the standard theoretical tool applied to the study of the electronic structure 
of metal surfaces has been Density- Functional Theory (DFT). 4 As in the original work of 
Lang and Kohn, most of the subsequent investigations have applied the Local-Density Ap- 
proximation (LDA) of DFT, or some of its semi-local variants (GGA, meta-GGA, etc.). 
This approach has been highly successful, and routinely yields good results for global sur- 
face properties such as work functions, surface energies, crystal-structure relaxation and 
reconstruction, etc- 

At a more basic level, however, some problems still remain to be solved, concerning for 
instance the asymptotic behavior of the exchange-correlation (xc) potential of the widely 
used Kohn-Sham (KS) approach to DFT. In the LDA, this potential decays exponentially 
when evaluated in the vacuum region, instead of the expected image-like oc — e 2 / z behavior. 6 
This qualitative failure of the LDA xc potential translates to a similar failure of the position- 
dependent xc energy per particle, e xc (r), which is defined through 7 



with -E xc [n] being the xc-energy contribution to the universal energy functional of DFT, 
and n(r) representing the electron density. Three aspects of Eq. ([1]) are worth emphasizing: 
(i) it represents the basic expression for the LDA, in which the exact e xc (r) of an arbitrary 
inhomogeneous electron system is replaced at each point r by that of a homogeneous electron 
gas at the local density n(r), and for a plethora of generalizations of the LDA,- (ii) since 
^(.[n] can be split as the sum of exchange (^[n]) and correlation (_E c [n]) contributions, 
one can write s xc (r) = e x (r) +e c (r), and (in) the position-dependent xc energy per particle 
e xc (r) entering Eq. (1) is not unique. One can always add to e xc (r) an arbitrary function 
e(r) with the condition that weighted by the electron density n(r) integrates to zero. Here 
we have chosen e xc (r) to represent the interaction between a given electron at r and its xc 
hole. 




(1) 
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The goal of this work is to provide exact analytical and numerical calculations of the 
exact-exchange e z (r) for jellium slabs and the semi-infinite (SI) jellium. In particular, we 
analyze the asymptotic behavior of the exact e x (r) in the vacuum region far away from 
the surface, and we find that there is a qualitative difference between e^ lah (z — > oo) and 
£3(2 — > 00): both exhibit an image-like behavior of the form —ae 2 /z(a > 0), but with a 
coefficient a that while in the case of jellium slabs is universal and equal to 1/2 in the case of 
a semi-infinite jellium depends on the density of the bulk material and only approaches 1/2 
in the low-density limit. The results reported here should help to settle the still controversial 
issue of the asymptotic behavior of the position-dependent xc energy per particle and KS xc 
potential at metal surf aces. 

Besides, being the results presented here exact at the exchange level, they should also 
serve as a benchmark against to which DFT xc calculations could be confronted, and hope- 
fully improved, once reduced to their exchange-only version. In this context, very recently 
Luo et alM have used a HF scheme to report self-consistent calculations of the surface energy 
and work function of jellium slabs. 

The rest of the paper is organized as follows: In Section II we present the general theo- 
retical background for both jellium slabs and the semi-infinite jellium, and we derive exact 
analytical expressions for the position-dependent exchange energy per particle in the vac- 
uum region far away from the surface. Numerical calculations that are valid at all positions, 
from the bulk region to the vacuum, are reported in Section III. Section IV is devoted to 
the conclusions. 

II. JELLIUM SLABS AND THE SEMI-INFINITE JELLIUM: THE EXACT 
ASYMPTOTIC BEHAVIOR 

A. Jellium slabs 

In the case of jellium slabs, with the discrete character of the positive ions inside the slab 
being replaced by a uniform distribution of positive charge (the jellium background), the 
positive jellium density is 

nf ah (z) =n9(-z)9(d + z), (2) 
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which describes a slab of width d, number density n,— and jellium edges at z = —d and 
z — 0. 9(z) represents the Heaviside step function: 9(z) — 1 if z > and = if z < 0. 
The size of the slab is infinite in the x — y plane. The jellium slab is taken to be invariant 
under translations in the x — y plane, so the KS eigenf unctions can be rigorously factorized 
as follows 

¥>i,k(r) = — r=£i(z), (3) 



where p and k are the in-plane coordinate and wave- vector, respectively, and A represents a 
normalization area in the x—y plane. £i(z) are normalized spin-degenerate eigenfunctions for 
electrons in slab discrete levels (SDL) % [i = 1, 2, ...) with energies Ef, they are the solutions 
of the effective one-dimensional KS equation 



h 2 d 2 



+ V KS (z) - e, 



U*) = 0, (4) 



2m e dz 2 

with m e being the electron mass. It is important to remark here that the factorization of the 
3D wave-function as proposed in Eq. ([3]) is only valid for the case of a local potential, as is the 
case of the KS implementation of DFT. On the other hand, in the HF approximation the non- 
locality of the Fock potential introduces a coupling between k and i quantum numbers^ 2 . As 
a consequence, HF numerical calculations are more time-consuming than the ones presented 
here, either LDA, KLI— , or OEP.i^ 

The local KS potential Vyis(z) entering Eq. (j3J) is the sum of two distinct contributions: 

V KS (z) = V H (z)+V xc (z), (5) 

where Vg.(z) is the classical (electrostatic) Hartree potential, given by^ 

/oo 
dz' \z - z'\ [n slah (z') - n+(z')] . (6) 
-oo 

Here, n slah (z) is the electron number density 17 

occ. 

- slab w = iE(^) 2 i^wi 2 ' (7) 



where k % F = \j2m e (eF — £i)/h, and ep = ep{n, d) is the Fermi energy or chemical potential, 
which in turn is determined from the neutrality condition for the whole system by the 
condition Y1T° (^p) 2 = 2ir dn. V xc (z) is the nonclassical xc potential, which is obtained as 
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the functional derivative of the ax-energy functional E xc [n(z)]-M 

1 6E xc [n(z)] 

Vxc[z) = A 5n(z) ■ (8) 

Both for the slab and semi- infinite [d — > oo limit of Eq. (|2j))] geometries, and as a 
consequence of the translational symmetry in the x — y plane, Eq. (JT]) simplifies to 

oo 

E xc [n}=A J dzn(z)e xc (z), (9) 

— oo 

where e xc (z) is the position-dependent xc energy per particle at plane z. In the case of 
jellium slabs, the exchange-only contribution to Kj C [n] (which is originated in the Pauli ex- 
change hole with all other correlation effects excluded) is known to be given by the following 
expression:— 



oo oo 



Ef ab [n] = -2e 2 A^44 dz dz' Vt (z, z') Vj (z' , z)F tj (z, z'), (10) 
where <Pi(z,z') = and 



-oo — oo 



oo 



A ' ) 4vry p + ( z _ Z 'f ' 1 ; 





with J\(x) being the cylindrical Bessel function of first order. 20 

Comparison of Eq. (Q and (|T0|) yields the following expression for the exchange-only 
contribution to e^ b (z)\ 



£ Slab, 



oo 

z ) = - ^b(zj Z2 k F k F J dz'<p i (z,z')p j (z',z)F ij (z,z'), (12) 



which can be interpreted as the energy due to the interaction of an electron at z and its 
exchange-only Pauli hole. In order to demonstrate that the exact-exchange energy per 
particle e x lab (z) of Eq. (I12p represents indeed the interaction between an electron at z and 
its exact-exchange hole, we appeal to the following expression for the exchange-hole for our 
slab geometry 2 ^- 

. occ. 

h x (r; r + R) = — — — - £ fcj4 J x (p^) J 1 (p4)&(z+Z)*& , (13) 
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which represents the density of the exchange hole at point r + R (observational point) due 
to the presence of an electron located at r. Owing to the translational symmetry in the 
x — y plane, without loss of generality we have choose r = (0, z), r + R = (p, z + Z). Using 
Eq. f[T5j) . and defining z' = z + Z, Eq. (fT2|) may be rewritten as 

eT\z) = $ ! dp I dz> M== > (14) 



z - z'V 



which justifies the physical interpretation of e^. lab (z) as the interaction energy of an electron 
located at z and the "charge distribution" given by h x (z;p,z'). Equations (IT5|) and (jHJ) 
can also be used as a sort of alternative definition of the e^. lah (z) investigated in this work, 
as they solve the non-uniqueness of (z) which results from its definition through the 
exchange-only version of Eq. Q.— ^ 



1. Single occupied slab discrete level 

In order to obtain the asymptotic behavior of e^} ab (z) at z — > oo, we first restrict our 
analysis to the case where there is one single occupied SDL.— In this special case, in which 
i — j — Eq. ((12} yields [see also Eqs. © and (fTTD ]: 

oo oo 2 

e ?f(*) = - e 2 / ^' |6(^')| 2 f ^ , (15) 

XA J J P y/ P 2 +(z-z') 2 

— oo 

or, equivalently (see Appendix): 

— oo 

with Ii and L\ being the modified Bessel and Struve functions, respectively.— 

We note that Eq. ffl6|) is valid for all z, both inside and outside the jellium slab. Also, 

the cancellation of n slab (z) which occurs in passing from Eq. (fl2l) to Eq. (JTB1) allows for the 

numerical calculation of e^ l 'f°(z) for arbitrarily large values of z. 

a. Asymptotic behavior. For the slab geometry, it is permissible (and rigorous) to take 

the asymptotic limit kp \z — z'\ ~ z kp ^> 1, although the integral over z' runs from — oo 

to +oo. This is due to the fact that for a given z the main contribution to the integral in 

Eq. (fTBj) comes from values of z' inside the slab (—d < z' < 0), as £%(z) decays exponentially 

one or two A^'s from each jellium edge. 

7 



h{2k F \z- z'\) L^kplz- z'\) 



kp \z 



k\ \z 



(16) 



At this point, we define F(x) = l—Ii(2x)/x+Li(2x)/x and use the asymptotic expansions 
of Ji and L\ in the limit x ^> 1. We obtain 20 

^>D-i-s + ^-- ' (17) 

Hence, as zk F ^> 1, Eq. ( fi~6l) yields the following asymptotic behavior: 

>— s( 1+ t + ? + -)< < 18 > 

where r s ) = z 1 ^, r s ) — 2/ [7rfc|i(d, r s )] and 71 (d, r s ) = z 2 (d, r s ) — 4 z 1 / [nk F (d, r s )], mean 
values being defined here as O* = J £i(z)*0(z) &(z) dz. The r s dependence of ~z l (d,r s ) and 
z 2l (d,r s ) comes from the self-consistent KS wave-functions £i(z), which for a given d are 
different for different values of the slab density dictated by r s . For details on the derivation 
of Eq. (jTBl) from Eq. ( fl~6l) . we refer to the Appendix. 

2. General situation 

For the general situation where more than one SDL is occupied, we obtain the asymptotic 
limit of Eq. (112p by using the fact that for z — > 00 (i) the electron density is dominated by 
the slowest decaying KS orbital, which corresponds to the highest occupied SDL (z = to), 
and (ii) the numerator of Eq. (fi~2l) is dominated by the term i = j = m, since all £j(V) with 
i m decay exponentially two or three A^'s from each jellium edge. Hence, in the vacuum 
region far away from the surface we find: 

n s ^(z^ oo) ^^l\U(z)\ 2 (19) 



2tt 



and 



00 

-,2 



47re f 

ef ah (z -> 00) -> - 2 / d2Vm(2, J)fmtf, z)F mm (z, z'), (20) 

or, equivalently [see Eq. ffTTl) ]: 



— oo 



A| 2 [Jl(p^)] 



2 



b (^ - oo) - - e 2 / dz' |e m (z')| Z / ^ 77"" . (21) 

^ VP + l- 2 _ 2 ) 







Finally, following the same procedure as in the case of a single occupied SDL, we find: 
4 lab (^oo)--^(l + ^ + ^ + ...), (22) 



where f3 m (d,r s ) = z m (d,r s ) - 2/ [irk%(d, r s )) and j m (d, r s ) = z 2 (d,r s ) - 
iz m (d,r s )/[7rk^(d,r s )).^ 

This result represents a straightforward generalization of the result presented above 
[Eq. ([18]) ] for the case of one single occupied SDL. At this point, it is interesting to note 
that the leading contribution to e^ lab (z — > oo) — > — e 2 /2z can be easily obtained directly 
from Eqs. f[T5]) or ([2"T]) . by approximating the argument inside the square root by z (in 
the large z limit), and using the normalization of the KS orbitals £i(z) and the identity 

oo 

J dxJ x {xf/x = 1/2. 



By considering a slab of thickness sufficiently large to make the energy spectrum continu- 
ous, Solamatin and Sahni^ reached the conclusion that far away from the slab the so-called 
Slater potential Vs(z) [which is twice the exchange energy per particle: Vs(z) = 2e x (z)] de- 
cays as —e 2 /z 2 , in contrast with the asymptotic structure dictated by Eq. (]22l . This result 
is, however, not correct due to the fact that for a finite jellium slab (no matter how thick it 
is) the slab intrinsic discrete spectrum [corresponding to the eigenvalues £j entering Eq. (jlj)] 
can never be replaced by a continuous one. 

Equation ([22]) leads us to the conclusion that in the vacuum region of a finite jellium 
slab and at distances from the surface that are large compared to l/h™ (which is typically 
larger than the slab thickness d), e^ lah (z — > oo) — > — e 2 /2z, which is exactly half the 
corresponding KS exact-exchange potential V x (z — > oo) — > — e 2 /z.— Hence, as in the case 
of finite systems,— the Slater potential Vg(z) of jellium slabs [or, equivalently, twice the 
exchange-energy per particle e x (z)] embodies the asymptotics of the KS exchange potential 
V x (z). 

In contrast, Solamatin and Sahni^^ concluded that in the case of a semi-infinite jellium 
only half the Slater potential embodies the asymptotics of the KS exchange potential, i.e., 
V x (z — ► oo) = e x (z — > oo); but Nastos^i claimed that V x (z — » oo) = 2e x (z — > oo), so there 
is still something remaining to be clarified on this issue. Work along these lines is now in 
progress. ~ 

B. Two-dimensional electron gas 

The exchange energy of a strict two-dimensional (2D) electron gas can be obtained from 
that of a jellium slab with a single occupied SDL [Eq. (flOj) with % — j = 1], by first performing 
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the one-dimensional non-uniform scaling 29 

nf ab (z) = Xn slah (Xz) , (23) 

and then taking the limit as A — > oo. The scaling above preserves the total number of 
electrons. Noting that for a single occupied SDL the jellium-slab exchange energy takes the 
following form 

— oo — oo 

where 



we find: 



n^(z) = ^Uz)\\ (25) 

271 

Ef ee lim E^ h [nf b ] = -iV A e 2 4 ; (26 ) 

A^oo cS7T 

where iV = A (A;]?) 2 /(27r) represents the total number of electrons. Previously, this scaling 
limit had been formulated in a different way, resulting in the much generous constraint that 
the exchange energy per particle in the 2D (A — > oo) limit should be greater than — oo.— 
It is interesting to note that the exchange energy functional as given by Eq. (|24p is an explicit 
functional of the density, which is only possible in this single occupied SDL case, due to the 
simple (invertible) relation between density and wave- function, as given by Eq. ( |25l) . In 
the general, many SDL occupied case, Eq. (1231) is replaced by Eq. ([7]), the direct inversion 
from wave-functions to density is not feasible anymore, and the exchange energy functional 
is an explicit functional of the KS orbitals, but an implicit functional of the density, as in 
Eq. (HO]). 

We note at this point that the exchange energy of a strict 2D electron gas can also be 
obtained directly from Eq. (115p through the replacement ^i(z') — > a/ S(z'), with S(z') being 
the Dirac delta function, and taking z = 0: 

oo 

Ef = -Ne* \% [J x (p4)] 2 = ~N± e 2 ^. (27) 







Either from Eq. ( |26l) or ( |27l) . we find for the exchange energy per particle of the strict 2D 
homogeneous electron gas the well-known result £~ 2D = _E 2D /N = — (4/37i)e 2 kp .— 
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C. Semi-infinite jellium 



In the case of a semi-infinite jellium, a half-space filled with a uniform distribution of 
positive charge (the jellium background), the jellium density is 

nf(z)=n6(-z), (28) 

with the jellium edge at z = defining the surface of a metal. As in the case of jellium 
slabs, the semi-infinite jellium is invariant under translations in the x — y plane, so the KS 
eigenfunctions can be factorized as follows 

— < r >=Sf' (29 » 

where p and k are the in-plane coordinate and wave-vector, respectively, and A (L) repre- 
sents a normalization area (length). ^k z ( z ) are spin-degenerate eigenfunctions for electrons 
with a continuous energy spectrum e kz = Vks(~ oo) + (hk z ) 2 /2m e (k z is a continuum quantum 
number). They are the solutions of the effective one- dimensional KS equation 

h 2 d 2 



4s 



+ Vks (z) - e kz 



&,(*) = 0. (30) 



2m e dz 2 

The KS potential Vks (2) is given by Eq. (j3J), as in the case of a jellium slab but with the 
slab electron density of Eq. ([7]) being replaced by the SI electron density 

kp 

= / (k 2 F - k 2 z ) \£ kz (z)\ 2 dh. (31) 



-k 



F 



In the case of a semi-infinite jellium, the position-dependent exchange energy per particle 
at plane z is given by the following expression: 

kp kp 00 

£ *( z ) = - 27r S(z) I dh J <(4-^ 2 ) 1/2 (4-A?) 1/2 / dz'^(z,z')^ z (z',z)F hX (z,z'), 

—kp —kp —00 

(32) 

where ^kX z i z ') = ik-A z )*ik-A z ')i anc l ^k z k' K z i z ') ls °f ^ ne f° rm of Eq. (fiTi) but with k l F (kp) 
being replaced by [kp — k 2 ] 1 / 2 ([kp — k' 2 } 1 / 2 ). 

a. Asymptotic behavior. The derivation of the asymptotic limit of e^(z) is more del- 
icate than in the case of jellium slabs, due to the fact that in the present case we have 
a continuous energy spectra. Hence, the crucial argument that we have used to derive 
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the asymptotic behavior for jellium slabs, concerning the fact that in that case the high- 
est occupied SDL dominates in the vacuum region far away from the surface, is not so 
transparent when the spectrum is continuous. Stated in other words, contributions to the 
asymptotic position-dependent exchange energy of Eq. (1521) come indeed from values of k z 
and k z that approach kp, but not necessarily only from the highest occupied value, i.e., from 
k z — — k z — — kp. 

The asymptotic behavior of Eq. (1321) was analyzed by Nastos^- by assuming that at 
z — > oo the KS potential Vk${z) takes the image-like form Vks(<£ —* oo) — > — «ks e2 / z , with 
«ks positive, but otherwise arbitrary. One finds that in the vacuum region far away from 
the surface [z — > oo) the KS orbitals ^ can be expanded with respect to the KS orbital at 
k z — kp as follows^ 1 ^ 

Uz - oo) - i kF {z - oo)e- Q2 ^- fe ) , (33) 

with 



i kF {z -> oo) oc e ~W 2m e w / h \2zVzW) aKs/ ^ 2mea2 ° w/h2 , (34) 

a standing for the square root of the ratio between the Fermi energy and the work function 
W {a 2 = e F /W). By introducing Eqs. ([33]) and (JMD into Eqs. (f3"T|) and (1521). one finds±i^ 

l«^°°>l'- (35) 

and 

Furthermore, the asymptote of Eq. ( |36i) does not depend on the actual form of ^k F {z — > oo). 
This is due to a cancellation, when z is large, of the orbitals £a; f (z — > oo) entering the 
numerator and denominator of Eq. (|32|) . This is the reason why the leading term in the 
expansion of £%(z — > oo) is independent of axs- That means that the result remains valid 
also in the absence of this image-like contribution, i.e., even assuming that the KS potential 
Vks(^) entering Eq. (Tj0) decays exponentially as z — > oo. 

The asymptotic behavior dictated by Eq. ( |36l) was obtained independently by Solamatin 
and Sahni using a somehow less general, but otherwise quite different approach .-^^ They 
approximated the KS potential VksW entering Eq. fl30l) by a finite-linear-potential model 
at the interface region, and used the corresponding orbitals in each of the three regions 
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where the model was defined: trigonometric functions (in the bulk region), Airy functions 
(near the surface), and exponential decaying functions (in the vacuum). As the only thing 
that matters in obtaining the asymptote of Eq. (I36p is the correct expansion of the KS 
orbitals with respect to the KS orbital at k z = kp [as given by Eq. f[3"3"j) ]. and with this 



genera 
and Id 



expansion being fulfilled also within the finite-linear-model potential used in Refs. [9 
Solamatin and Sahni obtained Eq. (I36I) which is valid in general. It is also worth 
of address the fact that the result of Eq. fl36l) is in contrast with the asymptotic behavior 
e x (z — > oo) — > — e 2 /Az that one obtains for the exchange energy per particle in the case 
of the Airy edge electron gas.— This is due to the fact that the asymptotic behavior of the 
solutions of the Airy edge gas is different from the one given by Eq. (13"3"j) . as for this model 
the potential increases linearly with distance in the vacuum region, instead of approaching 
a constant value. An analysis of the so-called Pauli and lowest-order correlation-kinetic 



components of the exchange energy per particle e x (z — > oo) can be found in Ref. 
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III. NUMERICAL RESULTS 

All the numerical calculations presented below have been carried out by ignoring all 
correlation effects (beyond the Pauli exchange). By this we mean that the xc potential V xc (z) 
entering Eq. (J3J) has been replaced by the exchange-only contribution V x (z), disregarding 
V c (z), both for jellium slabs and for the semi- infinite jellium. In the case of jellium slabs, 
V x (z) and the corresponding KS orbitals of Eq. (@J can be obtained through the solution of 
the discrete version of the x-only optimized effective potential (OEP) method, as given for 



example by Eqs. (14) or (20) of Ref. |36|. This code is feasible and we have at our disposal these 
self-consistent exact-exchange (OEP) KS orbitals and exchange potentials.—^ Alternatively, 
an approximate way to obtain V x (z) and the corresponding KS orbitals (the exchange-only 
LDA orbitals) is to replace the actual exchange potential V x (z) by the exchange potential 
of a uniform electron gas at the local electron density n(z), i.e., V x DA (z) = — [Qn(z)/n] 1 ^ 3 
[hartrees] . 

Fig. 1 shows a comparison of (i) the well-known LDA exchange energy per particle 

1 /3 

£x,lda{z) = — (3/47r) [37r 2 n(z)] [hartrees], with n(z) being the self-consistent electron 
density obtained with the use of either exchange-only LDA orbitals (wide solid line) or 
exact-exchange (OEP) orbitals (wide dotted line), with (ii) the exact-exchange energy per 
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FIG. 1: ef ab (z) for r s = 2.07 and slab width d = 4.3 X F . The curves denoted "Slab" have been 
evaluated from Eq. (|12p . by using exchange-only LDA orbitals (solid line) and exact-exchange 
(OEP) orbitals (dotted line). The curves denoted "LDA" have been evaluated from the well-known 
LDA formula e^LDA^) = — (3/47r)[37r 2 n(z)] 1//3 [hartrees], with n(z) being the exchange-only self- 
consistent electron density obtained with the use of exchange-only LDA orbitals (wide solid line) 
and exact-exchange (OEP) orbitals (wide dotted line). Inset: enlarged view of the bulk region near 
the surface. The bulk value of e x for this electron density is — 0.22134 (e 2 /ao). 

particle e^ lab (z) of Eq. (I12p obtained by using, as before, either exchange-only LDA orbitals 
(solid line) or exact-exchange (OEP) orbitals (dotted line). It is important to note that 
while both alternative evaluations of ^.lda^) fail badly in the vacuum region where the 
actual exchange energy per particle exhibits an image-like asymptotic behavior, the use of 
LDA orbitals in Eq. ffl2l) results in an exchange energy per particle (solid line) that on the 
scale of the figure is nearly identical to the exact fully-self-consistent result (dotted line). 
Small differences introduced by the use of LDA orbitals (see the inset) are mostly localized 
in the bulk region near the surface, where Friedel-like oscillations appear to be too weak in 
this approximation. 
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The numerical methods that we have used to obtain exact-exchange (OEP) orbitals and 
exchange-only LDA orbitals suffer from instabilities in the vacuum region far from the sur- 
face. As in this region exchange-only LDA orbitals are stabler than their OEP counterparts 
and the results presented in this work do not depend significantly on whether exact-exchange 
(OEP) or exchange-only LDA orbitals are used in Eq. (fT2j) (see Fig. 1), all the calculations 
presented below have been obtained with the use of exchange-only LDA orbitals. We em- 
phasize, however, that this is not a crucial approximation, and that the magnitude of the 
error that it introduces is given by the almost indistinguishable difference between the full 
and dotted lines in Fig. 1. 

The numerical self-consistent calculations presented in Fig. 1 and in the remaining of this 
section (which are all obtained from either Eq. (|P2|) or Eq. (I32|) with the use of exchange- 
only LDA orbitals) have been performed as follows. For jellium slabs, two infinite barriers 
have been located in the vacuum region far enough away from the two surfaces, in such a 
way that all the numerical results be independent of their precise location, 37 and the KS 
equations have been solved through a straightforward discretization in real space along the 
one-dimensional coordinate z. In the case of the semi-infinite jellium, the KS equations have 
been solved by following the general procedure introduced by Lang and Kohn.- This consists 
of defining three regions for the solution of the KS equation: far-left (bulk region), central 
(a few Ai?'s to the left and to the right of the jellium edge), and far-right (vacuum region). 
In the bulk region, the KS eigenfunctions are taken to be of the form £k(z) = sin(fcz — 7^), 
where 7^ are phase shifts, and this fixes an overall normalization constant. In the central 
region, we define a mesh of N points between z\ and z^ {z\ < zn), the first point z\ being 
chosen far enough from the jellium edge in the bulk so that the Friedel oscillations can be 
neglected, and the outer point zn being chosen to be far enough from the jellium edge into 
the vacuum so that the effective one-electron potential is negligibly small. Since Vks(^) ~ 
for z > zn, the orbitals can be approximated as = ae~ k * z where a is a constant 

and k* = (— 2m e €k z /h 2 ) 1 ^ 2 . The KS orbitals at the mesh points are calculated by using the 
Numerov integration procedure.— As in the vacuum region the orbitals follow exponential 
form, it is numerically most stable to integrate them inwards, so the Numerov integration 
procedure in this case is given by: 



2 + 10h(V KS (z i )-e kz ) 
1 - h(V K s(zi-i) - £k z ) 



ik{zi) 



1 - h(V K s(z i+1 ) - e fc J 
1 - h(V K s{zi-i) - £&,) 



(37) 
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FIG. 2: Position-dependent exchange-energy per particle for a thin jellium slab with r s = 2.07 
and d = 0.3 Xf ~ 2ao- Solid line: full numerical calculation of Eq. (|12p . Dashed-dotted lines: 
asymptote of Eq. (fT8|) ; the dashed (dotted) lines represent the asymptote of Eq. ([18]) with the last 
term (last two terms) neglected. Inset: enlarged view of the asymptotic region. 

where h = (z«+i — -2j) 2 /12. Finally, matching the KS orbitals in the central region with the 
corresponding analytical expression in the bulk region [£k( z ) = sin(kz — 7^)] determines the 
constant a. 

A. Jellium slabs 

In Fig. 2, we consider a thin jellium slab with r s = 2.07 (corresponding to the average 
electron density of Al) and d = 0.3 A^. This slab contains one single occupied SDL, so that 
we compare our full numerical calculation of Eq. ( fT2|) (solid line) with the asymptote of 
Eq. ( ITS!) (dashed-dotted line). We see that e^ lab (z) reaches Eq. (TTg|) at about one from 
the jellium edge and reaches the asymptote —e 2 /2z at a few Fermi wavelengths (~ 5 — 6Xf) 
from the surface. 
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FIG. 3: Same as Fig. 2, but for a slab with d = 4 \p, with nine SDL occupied. Full thick 
line, e^} ah (z) from Eq. (|12|) ; dotted, dashed, and dashed-dotted lines, asymptotic expansions from 
Eq. (|22p . The contribution from each pair (i,j) of occupied SDL to the total e x slab (z) are repre- 
sented with thin full lines, except for the last contribution (i = j = 9). Inset: enlarged view of the 
asymptotic region. 

In Fig. 3, we consider a jellium slab with r s = 2.07 and d = 4Xp. For this particular case, 
nine SDL's are occupied, i.e., eg < Sf < £io, so we compare our full numerical calculation of 
Eq. ({121) (solid line) with the asymptote of Eq. (|22|) (dashed-dotted line). The main message 
of this figure is that (i) m the vacuum region far away from the surface e^ lah (z) is dominated 
by the term i = j = m (dashed-dotted-dotted line), and (ii) e^ lab (z) reaches the asymptote 
—e 2 /2z only at a distance from the jellium edge of several Fermi wavelengths (~ 8\p)- We 
remind here that point (i) above was our main assumption in the derivation of the slab 
asymptotic limit of Eq. (!22j) . This assumption is fully justified after the numerical results 
shown in Fig. 3. 

At this point, with a few algebraic manipulations, we rewrite the asymptote of Eq. (|22|) 
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FIG. 4: Parameter zf ah (d, r s = 2.07, z) from Eq. ([35]), as a function of d, in several approximations. 
Full line, zf ah (d,r s = 2.07, z^ oo)+d/2 = -2/vr/c™(d, r s = 2.07); squares, Eq. dUJ; circles, fit to 
Eq. ([38]) . Upper panel, k™(d) as a function of d. 

in the physically motivated image-like form 



ef ah (z -> oo) 



where 



a 



«r b = 1/2, 



Slab 



[z-z®«*(d,r a ,z)]' 



(38) 



(39) 



and z| (d, r s , z), which represents the location of the so-called image plane, results in 

4" « r„ ,) = M d, r.) + + O (1) . (40) 

As 2 — >■ oo, z^ lah (d, r s , z) reaches a finite value, given by 

zf ah (d, r s , z -> oo) -> /3 m (d, r s ) = ^ m (rf, r s ) - 2/^ (d, r s ) = -d/2 - 2/irk^(d, r s ). (41) 

In Fig. 4, we plot a comparison of Eq. (HT]) (solid line) with the image-plane position that 
we obtain by fitting our full numerical calculation of Eq. f[T2]) with the image-like Eq. ([38]) 
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and a^ lab = 1/2 (empty circles). For this, we have used a fit region of width 2\p centered at 
6Xf from the jellium edge in the vacuum. Differences between Eq. (j4"Tl) (solid line) and our 
numerical estimate (empty circles), which in the case of very thin films are negligible, are 
entirely due to the fact that the fitting of the numerical calculation must be carried out in a 
vacuum region that extends very far away from the surface. Empty squares correspond to the 
result of Eq. (14*0]) . including the correction to the leading term, and taking z = 6Xf (which 
is the average value of the fit region indicated above). According to Eq. (|4"T1) . z^ lah (d) + d/2 
is inversely proportional to k^(d), which exhibits an oscillatory behavior as a function of the 
slab width d (see the solid line in the upper part of Fig. 4) going to zero every time a SDL 
becomes occupied. Hence, the location of the image plane becomes infinitely negative every 
time a SDL becomes occupied, which results in the strong finite-size oscillations shown in 
Fig. 4. 

B. Semi-infinite jellium 

Now we focus on a comparison between our full numerical jellium-slab and semi-infinite- 
jellium calculations of the position-dependent exchange energy per particle (see Fig. 5). In 
the bulk, as d increases the slab calculations converge with the semi-infinite calculation 
(see the inset at the upper part of Fig. 5), and both slab and semi-infinite calculations 
approach in the bulk region far away from the surface the exchange energy per particle of 
a three-dimensional (3D) homogeneous electron gas, e 3 ^ / (e 2 / Oq) = -(3/47r)(97r/4) 1 / 3 /r s « 
— 0.22134. In the vacuum, however, there is always a region far enough away from the 
surface where the jellium slab and the semi-infinite jellium behave differently: while all slab 
calculations converge to an image-like behavior of the form of Eq. fl38l) with a^ lab = 1/2, the 
semi-infinite e^(z —>■ oo) exhibits an image-like behavior in agreement with Eq. (1361) (see 
Fig. 6). Fig. 5 also shows that as the width d increases the slab e^ lah (z) coincides with the 
semi-infinite £^(z) in a wider vacuum region near the surface (see the lower inset of Fig. 5). 

Fig. 6 shows a comparison of our full numerical calculation of Eq. (1321) with the 
asymptote of Eq. fl36|) (solid and dashed lines, respectively) for r s = 2.07. In this case, 
a = \Jep/W = 2.048 (as obtained from our exchange-only LDA self-consistent calculation 
of the work function W) and Eq. (f3"6"j) yields £%(z — > oo) — > —0.18622 e 2 /z (dashed line), 
which is in contrast with the asymptote of Eq. (12"2"1) (dotted line) that holds in the case of 
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FIG. 5: ef ab (z) for r s = 2.07, slabs with d = 
semi-infinite case, from Eq. (|32p . Upper inset: 
asymptotic region. 



4,6,8, and 12 X F . Dotted line, ef-(z) for the 
enlarged view of the bulk region. Lower inset: 



jellium slabs. The same comparison has been done for other values of r s , and we have found 
that our full numerical calculation is always very close (as in Fig. 6) to the asymptote of 
Eq. (pjl. 

Finally, we display in Fig. 7 our exchange-only LDA self-consistent calculation of the 
coefficient a = \Jep/W in a wide range of electron densities. It is important to note 
that the corresponding coefficient [tt + 2aln(a)] / [27r(l + a 2 )] (see the inset to Fig. 7) en- 
tering Eq. (|36|) is close to 1/4 at metallic densities (r s = 2 — 6). Fig. 7 also shows that 
only at extremely low densities the coefficient a approaches zero, thereby the coefficient 
[tt + 2aln(a)] / [27r(l + a 2 )] of Eq. (|36l) approaching 1/2. Hence, the asymptotic limits of 
el lah (z) and e^ l (z) only coincide in the low-density limit (r s — > oo). 



20 




FIG. 6: Asymptotic behavior of (z) for the semi-infinite case, and comparison with Eq. f)36|) and 
the asymptote form of Eq. (I22p . 

IV. SUMMARY AND CONCLUSIONS 

We have presented a detailed analysis of the position-dependent exchange energy per 
particle e x (z) at jellium slabs and the semi-infinite jellium. 

For jellium slabs, we have found that in the vacuum region far away from the surface 
e^ lab (z — > oo) — > — e 2 /(2z), independent of the bulk electron density. This is the equivalent 
to the well-known result ^(r — > oo) — > — e 2 / (2r), which holds in the case of localized finite 
systems like atoms and molecules.- The equivalence between these results is however not 
straightforward, since slabs have an extended character in the x — y plane, being "localized" 
only along the ^-coordinate. In the vacuum side of the surface there is a region where e^ la,b (z) 
coincides with e^ l (z) and this region increases as d increases. The fitting of our numerical 
calculations of ef- ah (z) to a physically motivated image-like expression is feasible, but the 
resulting location of the image plane [z^ la,h (d,r s , z)] shows strong finite-size oscillations. In 
particular, we have shown analytically that z^. lah (d,r s , z — > oo) = —d/2 — 2/nkp(d,r s ), 
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FIG. 7: a(r s ) = ' EF(r s )/W{r s ) versus r s , for the semi-infinite case. Inset, the coefficient 
- [vr + 2aln(a)] / [2vr(l + a 2 )] of Eq. ([36]) versus r s . 

k™(d,r s ) being a signature of the energy of the highest occupied SDL with respect to the 
Fermi level. 

For a semi-infinite jellium, we have found that our numerical calculations agree well with 
the analytical asymptote [see Eq. (1361) ] obtained in Refs. and 
the slab asymptote — e 2 /2z only in the extreme low-density limit (r s 



33 . which approaches 
— > oo). 



We attribute the qualitatively different behavior of s% 



Slab/ 



oo) and ez l (z — > oo) to 



the fact that these asymptotes are approached in different ranges. While in the case of the 
semi-infinite jellium the asymptote is reached at distances z from the surface that are large 
compared to the Fermi wavelength (the only existing length scale in this model), for slabs 
the asymptote is reached at distances z from the surface that are large compared to 1/k™ 
(which is typically larger than the slab thickness d). For thick slabs with d» \p (Xp being 
the Fermi wavelength), e^ lah (z) first coincides with e^(z) [dictated by Eq. (36) at z » Xp\ 
in the vacuum region near the surface (see Fig. 5), but at distances from the surface that 
are large compared to 1/kp 1 , e^ lab (z) turns to the slab image-like behavior of the form of 
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Eq. (22) [or, equivalently, Eq. (38) with a^ lab = 1/2]; in the limit as d — > oo (i.e., when the 
jellium slab becomes semi- infinite), £ x h (z — > oo) coincides with £ x (z — > oo) everywhere. In 
the low-density limit, where A^ — > oo, the condition d >> A^ is never fulfilled and £^ lab (,2) 
reaches (at z » 1/k™) one single asymptote: the slab image-like behavior of the form of 
Eq. (22) [or, equivalently, Eq. (38) with a^ lab = 1/2], which turns out to coincide with the 
semi-infinite-jellium asymptotic behavior dictated by Eq. (36). 

Finally, we note that as e xc {z) = e x (z) + e c (z), the same conclusion is expected to be 
valid for the exchange contribution to the position-dependent xc energy per particle. Recent 
developments concerning the asymptotic behavior of the correlation contribution to the KS 
exchange correlation potential V xc (z) of a semi-infinite jellium can be found in Ref. 
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APPENDIX A: DERIVATION OF ASYMPTOTIC EXPRESSIONS 

Here we will shown in detail how to go through Eqs. f 1 1 5 j) - (TT5T) in the text. Starting from 
Eq. (ITo"]) . one first notice that^ 



dx Jf(ax) 



x 



^[x^Ttf %l 



■ _ h(2a\y\) Li(2 a |y|) 



a\y\ 



(Al) 



with 1\ and L\ being the modified Bessel and Struve functions, respectively. Substitution 
of Eq. (1 A II) in Eq. lIToT) yields at once Eq. (1T6l) . Now, in the asymptotic limit,— 



Lx(x > 1) -> h(x > 1] 



71 X 



2 • • • 7 



(A2) 
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which inserted in the definition of the function F(x) yields Eq. (fl7|) . Substitution of this 
expansion for F(x) in Eq. ffTB"]) leads to the expression, 



.Slab/ 



OO 



I z — z' I 



1 - 



7T fc^, | Z — Z' 



+ 



2{k\\z-z>\y 



— — I dz 



— oo 
AI2 



+ 



z — z 



dz 



ACi(z 



AI2 



Z — Z 



/|2 



2; — z 



/ 1 4 



(A3) 



(A4) 



Expanding the denominators of Eq. ( 1A4I) in the large z limit, the different contributions 
in Eq. (fl8l) arise. For instance, the leading contribution — e 2 / (2z) comes from the first 
term on the r.h.s. in Eq. (1A4|) . with \z — z'| approximated by z. It is interesting to note 
that one important feature of this result, as is its material and slab-size independence, is 
consequence (in this context) of the normalization of the SDL wave-functions. On more 
general grounds, and returning to the alternative definition of e^ lab (z) given in Eq. f fT4l) . this 
is more physically understood as a consequence of that the integral of h x (z; p,z + Z) over all 
possible "observational" coordinates (p, Z) is exactly —1.— The next term in the expansion, 
proportional to (3 1 and with decay z~ 2 , is obtained from the sub-leading contribution of the 
first term on the r.h.s. in Eq. (1A4j) . together with the leading contribution from the second 
term. To the order explicitly displayed in Eq. ( fl8|) . no contribution arises from the last 
(third) term in Eq. (1A4I) . as the leading contribution coming from this term to e^f^^z — ► 00) 
is of the order z~ 4 . While this analysis has been performed for the single-occupied SDL 
case, it also applies to the general case where more than a SDL is occupied, as explained in 
Section II.A.2. 
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